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A detailed description of the activation-relaxation technique (ART) is presented. This method 
defines events in the configurational energy landscape of disordered materials such as a-Si, glasses 
and polymers, in a two-step process: first, a configuration is activated from a local minimum to a 
nearby saddle-point; next, the configuration is relaxed to a new minimum; this allows for jumps 
over energy barriers much higher than what can be reached with standard techniques. Such events 
can serve as basic steps in equilibrium and kinetic Monte Carlo schemes. 
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I. INTRODUCTION 

Microscropic structural phenomena often proceed on 
time scales remarkably long compared to those of the 
atomistic oscillations. This is the case, for example, for 
glassy materials where microscopic dynamics takes place 
over time scales orders of magnitude larger than that as- 
sociated with the natural atomistic time scale, set by the 
phonon frequency of typically 10 13 Hz. Such a discrep- 
ancy is best understood from the configurational energy 
landscape: the system finds itself in a deep minimum sur- 
rounded by energy barriers which are many times larger 
than its temperature. Only rare fluctuations of thermal 
energies will allow the system to jump over a barrier 
and move to a new minimum. Typically, the rate for 
such jumps decreases exponentially with increasing bar- 
rier height, and may reach macroscopic values — of the 
order of seconds or more, rendering the study of these 
phenomena rather difficult. 

These long time scales are especially prohibitive for 
numerical studies. Traditional methods for the study of 
structural relaxation are of two kinds: molecular dynam- 
ics (MD) and Monte-Carlo (MC). MD is based on the 
direct integration of the equations of motion. In order to 
ensure the stability of the solution, the integration step 
cannot be larger than a fraction of a typical phonon vi- 
bration, i.e., somewhere between 1 and 10 fs. Depending 
on the number of atoms, the interaction potential, and 
the speed of the computer, the total number of steps can 
reach 10 4 to 10 7 , which translates into a time-scale on the 
order of nanoseconds; this is still far from the experimen- 
tal time-scale for structural relaxation of glassy materials. 
Because of the nature of MD, improvements beyond the 
linear level arc particularly difficult to achieve. Recently, 
a promising scheme involving a mixture of transition- 
state theory and MD has achieved a significant speed-up 
in the simulation of a model systemEJ; it is however too 



early to say how successful this scheme will be for generic 
problems. 

The inherent limitation to the degree of structural re- 
laxation achieved with MD does not apply a priori to 
MC schemes. Here, the speed of structural relaxation is 
mostly determined by the nature of the attempted moves. 
Until now, most algorithms have used moves defined in 
real space, involving the displacement of either one or a 
limited number of_atoms. Single-atom moves are rather 
efficient in liquidsa; however, they are not as successful 
in reproducing the collective nature of structural relax- 
ation associated with the slow dynamics of glassy and 
amorphous materials. Algorithms with more complex 
moves exist: the bond-switching algorithm of Wooten, 
Winer and WeaireJ3 for instance, succeeds in producing 
some of the best continuous random network models of 
amorphous semiconductors. Such algorithms are however 
problem-specific, and their dynamics generally unphysi- 
cal. 

In lattice models like the Ising model, it is often possi- 
ble to move from microscopic events — single spin-flips 
in the traditional Metropolis and heat-bath Monte-Carlo 
simulations - to collective events determining the behav- 
ior over longer times — flips of clusters of spins. Doing so 
can lead to a substantial improvement in the speed of al- 
gorithms, especially near the critical temperature where 
the correlation length and thus the cluster size diverges. 
The cluster algorithm of Swendsen and WangQ for exam- 
ple, can increase the computational performance of the 
simulation by many orders of magnitude compared to 
single-spin flip algorithms. 

In this paper, we give a detailed description of a 
recently proposed method which introduces a similar 
change of paradigm for continuum-based models: from 
the microscopic single-atom displacements to collective 
moves which form the basis of the activated processes 
in glassy and amorphous materials. This method, the 
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activation-relaxation technique (ART) , has already been 
applied with success to amorphous semiconductors and 
metallic glassesH~Q. With a similar algorithm, Doye and 
Wales have studied thq-potential energy surface of small 
Lennard-Jones clusters!! 

An event in ART is defined as a move from a local en- 
ergy minimum = (xf\ ■ ■ ■ , ) to another nearby 

minimum = (S^ , . . . , xty ) following a two-step 

process mimicking a physical activated processes: 

i) the activation during which a configuration is pushed 

from a local minimum to a nearby saddle-point; 

ii) the relaxation which brings the configuration from 

this saddle-point to a new local minimum. 

By defining the moves in the 3iV-dimensional space 
controlling the dynamics of relaxation - the configura- 
tional energy landscape - ART removes any constraint 
on the type of real space moves allowed. This is par- 
ticularly important in disordered and complex materials 
where events can involve very complex local or collective 
rearrangements that are hard to foresee. 

This paper is organized as follows: we first present the 
activation-relaxation technique. The following section 
discusses the implementation of the algorithm. We finally 
show examples of events in amorphous silicon (a-Si) and 
silica glas (g-Si02). 

II. THE ACTIVATION-RELAXATION 
TECHNIQUE 

In many materials and systems, the dynamics can be 
accurately described as a sequence of metastables states 
separated by energy barriers high compared to fceT, the 
typical energy scale at the atomic level. Such metastable 
configurations will remain essentially unchanged on a 
time long compared with the natural time-scale set by 
lattice vibrations, and can be well identified by the 
atomic positions at zero K, i.e., by the local minimum of 
the configurational energy landscape. Knowledge of the 
distribution and properties of these local minima is suffi- 
cient for determining the thermodynamical properties of 
the system. To understand the dynamical properties of 
these materials, however, a knowledge of the rates con- 
trolling the jumps from one minimum to another is also 
necessary. 

To a first approximation, the dynamics in these ma- 
terials is determined by the activation energy, i.e. the 
energy needed to bring a configuration from the local 
minimum to a nearby saddle-point. Because of the ex- 
ponential nature governing the energy fluctuations, any 
event following another trajectory, with by definition an 
energy higher than that at the saddle-point, will be much 
less probable and can be safely neglected.^ For the sim- 
plest characterization of the non-equilibrium properties 
or dynamics of a disorder material away from the glass 



transition it is therefore sufficient to map the continuous 
configurational energy landscape onto a network formed 
by minima connected, via trajectories going through first- 
order saddle-points JlJ The current ART method provides 
a local prescription for exploring this simplified space 
through hops from a local minimum to another (events). 

By defining the events directly in the configurational 
energy landscape, which, as we have seen, fully deter- 
mines the dynamical and equilibrium properties of a ma- 
terial, ART becomes much less sensitive to the details of 
the real space configuration. Doing so, it refrains from 
defining a priori the type of atomic rearrangements lead- 
ing to structural relaxation. In effect, it is the system it- 
self which determines the appropriate atomic processes, 
in much closer agreement with real processes. Such a 
change in paradigm, from real to configurational space, 
is particularly necessary for the study of glassy materials 
where an unambiguous description of real-space config- 
urations in terms of neighbor lists, coordination defects, 
etc./ is generally impossible to give. ART is a priori 
blind to the details of real space configurations; all ART 
needs is a local and continuous description of an energy- 
landscape; discontinuous energy landscapes, as, for in- 
stance, in discrete spin models, cannot be differentiated 
and thus forces are not defined. Any continuous interac- 
tion potential however, from Lennard-Jones to LDA, can 
in principle be used with ART. 

As mentioned in the introduction, the activation- 
relaxation technique consists of two parts: a path from 
a local energy minimum to a nearby saddle-point — the 
activation; and a trajectory from this point to a new 
minimum — the relaxation. 

The relaxation to an energy minimum poses no par- 
ticular challenge: it is a well-defined and well-behaved 
operation for which a number of efficient algorithms are 
available (see, for example, Ref. [Tl| ). 

The activation from a minimum to a saddle-point re- 
quires more care: to our knowledge, no theoretical frame- 
work exists that allows for finding the complete set of 
saddle-points around a local minimum. A number of 
works have been devoted to the study of finding the 
transition states in clusters and low-dimensional systems. 
Many of the techniques, however, start with the knowl- 
edge of both rrunirrjjum states and try to find the path 
connecting the two.E£l It is a very different problem to 
try to find a saddle-point with the knowledge of only 
one minimum. Most methods can be traced-, back to 
two techniques, the distuiauished coordinatell3o and the 
eigenvector-followingHOEa algorithms. Although these 
methods are generic, neither addresses the question of 
the generation of a complete set of saddle-points around 
a given minimum. 

In steepest-descent — or zero-temperature Langevin 
dynamics where the velocity is proportional to the force 
- all trajectories, including those starting at a saddle- 
point, lead to a local energy minimum. A naive approach 
to find the trajectory from a minimum to a nearby saddle- 
point would therefore be to retrace this path using a 
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time-reversed zero-temperature Langevin dynamics, or 
steepest- ascent algorithm. This fails, however, since us- 
ing steepest-ascent simply corresponds to inverting the 
sign of the total energy, in effect exchanging local minima 
with local maxima. Moreover, the minimum-energy tra- 
jectory leading from a local minimum to a saddle-point 
is an unstable trajectory for steepest-ascent; any pertur- 
bation sends the path away from the reversed steepest- 
descent trajectory. 

Within Newtonian mechanics a trajectory from a 
saddle-point to a minimum is also time-reversible: start- 
ing at a minimum with properly chosen velocities, one 
would be able to move up to any saddle-point. In con- 
trast to time-reversed Langevin dynamics, the trajectory 
cannot reach divergent parts of the configurational en- 
ergy landscape since the total energy is conserved. As 
with time-reversed zero-temperature Langevin dynam- 
ics, though, even a very tiny deviation can bring the 
system far away from the saddle-point. Sampling a 
very large number of initial random displacements and 
then targeting the least divergent trajectories, Dykman 
et al. could successfully find rthe saddle-points of a 
chaotic two-dimensional system.£3 If this approach can 
work for a simple energy function in low dimensions, 
such a hit-and-miss algorithm becomes hopeless in a large 
3N-dimensional space with a computationally expensive 
force to evaluate. 

At the saddle-point, all eigenvalues of the Hessian but 
one are positive. The energy landscape resembles a val- 
ley going down along the eigen directions correspond- 
ing to the negative eigenvalue. Leaving the saddle-point 
by steepest-descent we follow the floor of the valley to 
eventually arrive at a nearby minimum. This suggests 
immediately a local algorithm which should be more sta- 
ble than the steepest- ascent: to define a trajectory to a 
saddle-point, the configuration is moved in such a way as 
to minimize the force along all directions but the one cor- 
responding to the lowest eigenvalue. This eigenvalue is 
identified with the local bottom of the valley, and the con- 
figuration is moved against the force along this direction. 
A small displacement away from the bottom of the valley 
would be corrected for by the (3A^ — l)-dimensional min- 
imization, making the trajectory stable. Intuitively, this 
line and the path of steepest-descent should run mostly 
parallel; they are not identical though, and sometimes 
diverge. 

In most circumstances, this algorithm will converge to 
a saddle-point. Because we consider here the maximiza- 
tion along a single eigendirection, this algorithm will not 
lead to second- or higher-order saddle-points. This is in 
essence what was proposed by Cerjan and Miller for the 
location— of transition states in low-dimensional energy 
surfaces^, and what was used for an extensive study of 
a 13-atom LJ cluster by Doye and Walesll. 

Because of its N 3 requirements, this algorithm be- 
comes rapidly too computer intensive for realistic bulk 
systems, often demanding many hundreds of atoms with 
a costly energy function. We must therefore find another 



algorithm which does not require evaluation of the full 
Hessian matrix at each step. 

The current implementation of ART follows a modi- 
fied force vector G, obtained by inverting the component 
of the force parallel to the displacement from the cur- 
rent position to the local minimum r = X - M(°) while 
minimizing all other 3 AT — 1 directions: 



G = F — (1 + a)(F ■ f)f 



(1) 



where f is the normalized vector parallel to r, F is the 
total force on the configuration as calculated using an 
interaction potential, and a is a control parameter. This 
equation is applied iteratively until the force parallel to 
the displacement from the minimum F ■ f changes sign 
from negative to positive. Generally, the force perpendic- 
ular to the displacement decreases rapidly after a few it- 
erations, bringing the configuration close to the steepest- 
ascent trajectory. For a steepest-ascent path perfectly 
parallel to f , the modified force of Eq. ([!]) strictly sticks to 
the floor of the valley up to the saddle-point; for steepest- 
ascent trajectories perpendicular to r, the algorithm fails. 
From experience, such trajectories are rare and the algo- 
rithm generally converges to a saddle-point . 

Since moves are defined in the configurational energy 
landscape, vectors in Eq. (Q) laave 3N components both 
for the force and the positionjlla the displacement of the 
configuration from a local minimum to a nearby saddle- 
point may therefore involve any number of atoms — from 
one to all N atoms. 

In disordered networks, it is unlikely that the lowest 
eigenvalue of the Hessian matrix is degenerate. There are 
therefore always only two valleys stemming out of the lo- 
cal minimum, corresponding to the positive and negative 
direction of the lowest eigenvector. Thus, following val- 
leys from the minimum either along the lowest eigenvec- 
tor or the modified force leads to only two saddle-points, 
whereas a system typically has many, even thousands of 
saddle-points. Even worse, these two directions corre- 
spond, in bulk materials, to long-wavelength distorsions 
and do not lead to interesting events. Finding a way to 
avoid these directions can be a difficult task. 

One approach, taken by Doye and Wales for the study 
of a 13-atom Lennard-Jones cluster, is to select in turn 
each of the eigendirections of the Hessian .at the mini- 
mum and follow it to a nearby saddle-point a Since there 
are only 78 such directions, only a fraction of the many 
(« 10 3 , see Ref. Q saddle-points can be reached this 
way from the minimum; local information around the 
minimum is insufficient to locate all valleys leading to 
saddle-points. Moreover, the repeated calculation of the 
Hessian is an expensive operation for large systems. 

We propose a few approaches that do not require 
0(N 3 ) operations and work for a wide spectrum of cir- 
cumstances; these are discussed in section 
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Once a valley has been found, the situation becomes 
more straightforward, and we can use either of the algo- 
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rithms described above to follow the valley to the saddle- 
point. 

III. THE IMPLEMENTATION 

The implementation of the method poses no particu- 
lar conceptual or computational problems. The whole 
code, except for the force and total energy calculation, 
contains a few hundred lines at most. Its core consists 
of three parts: the escape from the harmonic basin, the 
convergence to the saddle-point, and the relaxation to a 
minimum. 



A. Escaping the harmonic basin 

The part of the algorithm that is most sensitive to 
details of the system studied is the escape from the har- 
monic basin; different approaches might have to be tried 
to find the most effective one. In general, open but stiff 
materials like amorphous semiconductors have a very 
small harmonic basin from which it is easy to escape. 
More compact materials — metallic glasses — or flop- 
pier ones — silica glasses — pose more problems. To 
ensure a proper sampling of events, any method for es- 
caping the harmonic basin that leaves out a significant 
fraction of the saddle-points should be avoided. 

The simplest way of escaping a harmonic basin is to 
make a random displacement away from the minimum, 
involving a single randomly chosen atom, a cluster, or all 
atoms. In our experience, for small systems they all lead 
to the same type of events; for larger systems, a global 
random displacement tends to induce many spatially sep- 
arated events which become difficult to disentangle. We 
therefore prefer a local displacement for systems of more 
than a few hundred atoms. 

A random direction generally has a sizable overlap with 
the softest elastic modes, and tends to fall back to these 
easily. We get better results by taking the escape direc- 
tion for the initial displacement along the force induced 
by a small random displacement; this procedure is es- 
sentially equivalent to applying the Hessian matrix to a 
random vector, resulting in a first-order suppression of 
the softest elastic modes. 

For small systems, where the Hessian can be obtained 
and diagonalized in a relatively short time, the softest 
modes can be removed directly from a random initial 
direction, or the initial displacement can be chosen along 
a linear combination of the stiffest eigendirections. This 
approach is rather computationally involved and cannot 
be reasonably carried out for systems with more than 100 
or 200 atoms. 

Once the initial direction is fixed, it is then followed 
until the passage of some threshold, indicating that the 
harmonic region has been left. This threshold has to be 
large enough to ensure that the trajectory does not fold 



back onto the softest direction while remaining insde the 
basin of attraction. We use a combination of two condi- 
tions for determining the point where the configuration 
has left the harmonic region surrounding the initial mini- 
mum: when the force component parallel to the displace- 
ment either stops increasing or when the ratio of this 
component to the perpendicular component is smaller 
then a given fraction, we consider the harmonic region 
left and the ART procedure as such begins. 

In the algorithm used in Refs. ||, || and |?|, no clear 
distinction was made between leaving the harmonic re- 
gion and convergence to the saddle-point; instead, an 
additional repulsive harmonic potential was introduced, 
which is added around the minimum with a strength A rep 
and a range r c : 

Erep = A rep (\r\ - r c f . (2) 

Although relatively efficient, this approach modifies the 
local energy landscape and introduces an artificial length 
scale r c in the problem. To reduce the impact of this 
additional length scale, one can re-initialize r c and A rep 
at random before each event. 

Currently we prefer to take as initial direction the force 
after a random displacement, and follow that direction 
until we leave the harmonic region, and then follow G as 
defined in Eq. (|f|) until the saddle-point is reached. 

B. Convergence to a saddle-point 

Convergence to the saddle-point cannot be achieved us- 
ing standard minimization techniques because the mod- 
ified force G as defined in Eq. ([!]) is not curl- free, i.e., 
it cannot be obtained from the gradient of a scalar. We 
therefore have to follow closely the direction of G until 
we reach the saddle-point, indicated by a change of sign 
in the component of the force parallel to f. Many sim- 
ple algorithms can readily be adapted for this purpose. 
Making small displacements in the direction of G is the 
most obvious choice for reaching a saddle-point. Such a 
crude method, however, is rather unstable and can easily 
enter into oscillations or severe slowing down. 

The conjugate-gradient (CG) algorithm provides an 
easy solution to this restriction by ensuring that the new 
displacement will be in a direction conjugate to the previ- 
ous onesJiil The line minimization along a direction h re- 
quired in the CG implementations of numerical packages, 
however, are based on the existence of a total energy - 
which cannot be defined. We replace it by a root-finding 
algorithm of G ■ h. In general, only a couple of force 
evaluations are necessary to reach that point. .— . 

The Levenberg-Marquardt (LM) aigorithma pro- 
poses a mixture of steepest-descent and a full-fledged 
second-order Hessian minimization technique. Away 
from the harmonic regime, the steepest-descent controls 
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FIG. 1. An event in the simulation of amorphous silicon. 
From left to right, the initial, saddle-point, and final config- 
urations are shown. The top and bottom rows correspond to 
different viewing angles of the same event. Dark atoms change 
their bonding environment during the event; light atoms are 
nearest-neighbours of the dark atoms. Activation energy: 
5.74 eV, energy difference from initial to final configuration: 
2.30 eV. 

the configuration reaches a stable energy. To obtain a 
low-energy configuration, we use the standard Metropo- 
lis algorithm where a new configuration is accepted 
with probability 1 if the energy is lower than that in 
the original configuration, otherwise with probability 
exp(— AE/hsT). The temperature as such is fictitious 
and we find that fcgT = 0.25 gives satisfactory-DDSults. 
As in Ref. |, we use a modified Stillinger-WebeiO inter- 
action potential with a three-body force twice the original 
value to remove the liquid-like features of the amorphous 
phase associated with the original SW. 

One event obtained in the relaxed structure is shown 
in Figure |l| from two difference angles. In the bottom 
representation, we can see how the configuration passes 
from three five-membered rings (initial) to one five- and 
one eight- mcmbered ring (final). In the process, four 
bonds are broken and four are created, preserving the 
total coordination, and the displacement incurred by the 
atoms is 2.3 A. This event has an activation energy of 5.74 
eV and the final configuration is 2.30 eV higher than the 
initial one. 

For silica glass, we use a 576-atom configuration re- 
laxed from the melt using molecular dynamicsE3. The 
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V. CONCLUSION 

By defining events directly in the configurational en- 
ergy landscape, the activation-relaxation technique pro- 
vides a generic approach to study relaxation in complex 
systems such as glassy and amorphous materials, poly- 
mers, and clusters. Real space moves are determined 
by the system itself and represent the most likely physi- 
cal trajectories followed during relaxation. ART is much 
less sensitive to the slowing down caused by increasing 
activation energy barriers than standard MC and MD 
approaches. 

Already ART has produced results which could not 
be achieved via other techniques: it. has produced 
well-relaxed samples of OrSiJa a-GaAs,EH3 Ni 80 P 2 oj3 and 
minimum-energy configurations of clusters of Lennard- 
Jones particlesfl The examples of events presented here 
demonstrate that ART can easily reach regions of the en- 
ergy landscape which are difficult to sample using more 
standard techniques. This paper provides the necessary 
description of the algorithm to allow for a rapid applica- 
tion of ART to a wide range of problems. 
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